Development  of  Algorithms  and  Prog^c'ams  for 
the  Ray  Radlotomography  of  the  Ionosphere. 


Final  (12-raonth)  report  (BOARD  Contract  SPG  93-4047) 


OORTEHT 


Introduction. 


19990211  031 


1 .  The  solution  of  the  direct  problem  of  radio  T/ave 

propagation  for  the  Ionospheric  Tomography.  8 

Description  of  the  program  for  calculation  the  phase  and  the  11 
phase-difference  for  arbitrary  positions  of 
the  receivers  and  the  satellite. 

2.  Design  of  the  different  versions  of  the  HT  operators  (matrices).  14 

Description  of  the  program  for  calculation  of  the  different  23 
versions  of  the  RT  matrices  (operators). 

3.  The  solution  of  the  Inverse  problems  of  the  phase  and 

phase-difference  RT.  25 

Description  of  the  program  for  solution  of  the  inverse  34 
problems  of  the  pliase  and  phase-difference  RT. 

4.  Analysis  of  the  influence  of  data  errors  and  noises  on  the  35 

reconstruction  results. 

References  s  ^2 


Figures  and  tables 


WETHIBUTION  STATEMENT  A 


Appiorcd  Iss  piab!!^ 
Dlstrlkiatiosj  UElJm!t®d 


Principal  investigator 


BUG  ©lii^IjXxY-'IKSlrEOTED  & 


Professor  V. Kunitsyn 

Reproduced  From 
Best  Available  Copy 


I 


1 .  Introduction. 

We  shall  consider  the  potentials  of  radiotomography  (RT) 
using  phase  and  phase-difference  measurements.  Task  of  the  ray  RT 
of  large-scale  structures  are  usually  formulated  as  follows: 
making  use  of  measurements  of  linear  Integrals  for  series  of  rays 
Intersecting  a  certain  area,  it  is  necessary  to  reconstruct  the 
structure  of  this  area.  Since  the  dimensions  of  large 
inhomogeneities  of  natural  (such  as  an  ionospheric  trough)  or 
artificial  origin  (traces  of  spacecraft,  technological  releases 
etc.)  are  of  hundreds  and  thousands  kilometers,  the  diffraction 
effects  in  the  case  of  VHP/UHP  probing  can  be  neglected  f1,2]. 

It  seems  to  be  not  possible  even  to  name  the  kinds  and  types 
of  emissions  and  waves  for  which  attempts,  at  least  theoretically, 
have  not  been  made  to  apply  tomographic  methods.  It  makes  no  sense 
here  to  list  all  of  the  variants  of  tomographic  approaches,  but 
perhaps  we  should  briefly  characterise  the  approaches  in  the 
closest  fields  of  tomography  of  geophysical  structures.  The 
methods  and  equipment  for  seismic  tomography,  where  the 
distribution  of  the  seismic  ’’slo-wness'*,  which  is  a  value  which  is 
the  inverse  of  the  wave  velocity,  is  reconstructed  from 
measurements  of  the  propagation  time  of  seismic  waves,  have  now 
been  developed  rather  well.  There  is  a  very  large  nimiber  of  works 
on  seismic  tomography;  we  will  only  cite  certain  surveys  [3-5]. 
Many  solution  methods  and  algorithms  have  been  proposed  in  the 
field  of  seismic  tomography  which  are  also  suitable  for  other 
types  of  waves.  Tomographic  studies  of  other  "spheres of  the 
Earth  are  also  carried  out.  Acoustic  tomography  of  the  ocean, 
where  the  propagation  time  of  acoustic  waves  is  also  measured,  la 
now  being  developed  actively  [6,71.  Works  on  radio  and  optical 
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tomograpliy  of  the  atmosphere  have  appeared  18,91. 

Proposals  on  variants  of  RT  of  the  Ionosphere  have  also  been 
offered.  A  tomographic  methods  for  determining  the  local 
attenuation  factor  as  a  function  of  geographic  coordinates  using 
data  on  the  integral  attenuation  factors  on  paths  covering  a  given 
region  with  a  sufficiently  dense  network  is  discussed  in  [10,111. 
iinltters  and  receivers  far  from  the  region  to  be  reconstructed 
provide  straight  propagation  paths.  It  was  suggested  to  use  the 
Radon  transformation,  which  is  practically  unfeasible  due  to  a 
small  number  of  rays,  besides,  the  question  of  the  dependence  of 
the  attenuation  factor  on  height  remains  to  be  settled.  A  variant 
of  RT  using  a  satellite,  where  reconstruction  of  the  electron 
concentration  distribution  from  the  TEC  on  a  series  of  rays  was 
proposed,  was  examined  in  [12,131.  The  data  acquisition  scheme  is 
shown  in  Pig.O,  where  three  receiving  stations  are  located  on  the 
plane  of  passage  of  the  satellite  and  the  rays  are  assumed  to  be 
linear.  The  schemes  for  tomography  according  to  linear  integrals 
for  various  types  of  waves  are  Indlstlguishable  from  one  another. 
The  approach  according  to  the  ray  integral  electron  concentration 
examined  in  [12,131  is  identical,  for  example,  to  seismic 
Interstitial  tomography  [141,  where  the  sources  of  seismic  waves 
are  placed  in  bore  holes,  while  the  receiver  is  moved  along  the 
surface.  More  recently,  methods  of  the  RT  of  the  ionosphere  using 
the  TEG  have  been  developed  by  many  investigators  [15-231. 

Prom  our  point  of  view  such  a  generalization  of  tomography 
using  linear  integrals  or  TEC  is  enough  obvious.  The  problem  is 
how  to  determine  this  absolute  phase  or  the  linear  integral.  Many 
Investigators  attempted  to  answer  this  question  and  discussed 
difficulties.  Here  a  nontrivial  point  is  that  when  determining  the 
absolute  phase  proport lanal  to  TEC,  the  problem  of  determining  the 


L 


3 


"Initial  phase"  arises  and  one  can  make  a  significant  error  and  be 
off  by  a  constant,  which  would  lead  to  inconsistent  tomographic 
data  and  render  the  reconstruction  unfeasible.  Unfortunately  the 
authors  of  earlier  publications  on  the  Ionosphere  tomography 
merely  presented  this  common  and  evident  Idea  and  failed  to  show 
the  possibility  of  the  reconstruction  in  the  presence  of  typical 
errors,  that  Is  they  did  not  model  the  Influence  of  errors  at  the 
initial  stages  on  the  reconstruction  results.  The  above  comments 
should  not  be  understood  as  diminishing  the  significance  of  the 
earlier  studies  on  Ionospheric  RT,  of  special  Importance  are 
studies  on  limitations  of  ionospheric  RT  and  resolution  limits 
Imposed  by  geometrical  and  sampling  limitations  in  Ionospheric  RT 
1 15, 163.  However,  the  authors  of  the  earlier  works  restricted 
themselves  to  simulating  the  possibility  of  reconstructing  the 
function  using  linear  Integrals,  which  was  of  little  use  since 
such  simulation  with  the  term  substitution  ("group  delays"  -♦  TEC, 
"seismic  slowness  distribution"  -*  "electron  concentration 
distribution"  and  so  on)  has  been  made  In  numerous  earlier  works 
on  seismic  and  other  kinds  of  tomography.  In  contrast  to  selsmlcs, 
however,  where  a  linear  Integral  (a  group  delay)  is  measured 
directly,  for  the  Ionosphere  there  are  techniques  for  only 
approximate  determination  of  linear  Integrals.  Here  a  specific 
error  appears  -  an  error  by  a  constant,  which  is  essential  for 
tomography.  In  direct  measurements  of  linear  Integrals  this  error 
has,  as  a  rule,  a  nolse-like  nature,  which  does  not  seriously 
affect  the  results.  As  discussed  earlier  £1,2,243  and  will  be 
sho?m  below  In  section  4,  the  phase  RT  leads  to  poor  results  In 
the  case  of  typical  errors  In  determining  the  absolute  phase  or 
TEC.  Therefore,  some  years  ago  £2,243  we  proposed  the 
ptiase-dlfference  RT  of  the  ionosphere. 


The  theoretical  foundation  of  the  ray  RT  la  the  well-known 
relations  [25,1]  for  phase  of  radlowaves  In  the  geometrical  optics 
approximation.  The  following  equality  determines  linear  Integrals 
from  the  electron  concentration  distribution  N: 

J  N  da  =  <|),  ^ 

where  X  Is  length  of  probing  waye,  r^  Is  the  classical  electron 
radius,  a)=kc,  k  Is  the  wave  number  In  the  free  space,  c  Is  the 
velocity  of  light,  /da  Is  the  symbol  of  Integrating  by  the  ray. 
Here  the  linear  Integrals  are  the  phase  difference  of  the 

field  being  measured  (E=Aexp(l<&) )  and  the  probing  wave  field 
(EQ=AQezp(l®^)).  Measurements  of  ^  are  realized  by  receiving 
signals  at  two  coherent  frequencies  from  navigation  satellites. 

We  Introduce  a  series  of  parameters  characterizing  the 
geometry  of  the  recording  scheme  in  polar  coordinates  (r,  a).  In 
Flg.O  (r^,  Oq)  -  are  the  coordinates  of  the  satellite,  (R,  a^)  - 
are  the  coordinates  of  one  of.  the  receivers  located  on  the  surface 
of  the  earth  (r=R);  p  is  the  elevation  of  the  satellite;  (|)  = 

(p-n;/2)  Is  the  angle  to  the  satellite  measured  from  vertical;  0  is 
the  center  of  the  earth;  axis  0-0’  Is  axis  of  the  polar  coordinate 
system.  On  the  basis  of  simple  geometrical  relations  for  any  point 
in  the  Ionosphere  with  coordinates  (r,  a)  located  at  the  distance 
of  1  from  the  receiver  the  foollowlng  equations  are  satisfied: 

^  r  _ R 

sln(a^-a)  sln(7i:/2+|3)  sln{7i:-(i3:/2+p)-(a^-<x))’ 

Prom  this,  we  obtain  an  equation  for  r(a)  of  the  straight 
(§=const)  ray 

r(a)  =  (R  cos  p)/(cos(p  +  -  a))  (2) 

and  the  relation  which  is  Inverse  of  It 


(3) 


R 

a(r)  =  +  p  -  arccos(-jr“  cos  p). 

The  relation  between  p  and  a,  r  follows  from  (2): 


tg  p  =  (cos(a^-a)  -  R/r)/(sln(a^-a)).  (4) 

Using  (3),  we  obtain  a  formula  for  an  element  of  the  ray 
of  length  da 

„  o  da  o  o  o 

do^  =  C1  +  r^( - ri  dr'^  =  - 5 - o"  ^  •  (5) 

dr  r  -  R  cos  p 

Then,  the  relation  for  the  measured  linear  Integral  (1 )  with 
respect  to  the  electron  concentration  will  have  the  foiro 

f  N(r,a)  r  dr 

I  N(r,a)  do  =  A,r_  —  - - .  (6) 

-*7  1^-  R^sln^q) 

In  the  place  of  the  polar  coordinates  (r,  a),  hereafter  It  is 
convenient  to  use  the  orthogonal  system  (h,  ’t):  h=(r-R)  Is  the 
height  above  the  earth’s  surface  and  t  =  oR  Is  the  ’’transverse” 
(horizontal)  distance  along  the  earth’s  surface  In  the  plane  of 
passage  of  the  satellite.  Here,  ray  equation  (2)  will  no  longer  be 
a  straight  line 

cos  p 

h(T)  =  R  [ - 11.  (T) 

cos(p+('t^-T)/R) 

Here,  (T:j,,h=0)  are  the  coordinates  of  the  receivers.  The  relation 
which  is  the  Inverse  of  (T)  Is  similar  to  (3). 

R 

t:  -  t:.  =  R  cp  -  arccos( —  cos  p)3.  (8) 

^  R+h 

In  this  case,  subject  to  (5),  the  linear  Integrals  of  type 


(1 )  have  the  form 


P(h,T)(R+h)(]li 

. =  KP.T.).  (9) 

Q  7  R^sln^p  +  2Rh  + 

Integration  with,  respect  to  the  connecting  receiver  i  (T^=a^R) 

to  the  satellite  is  replaced  according  to  (5)  hy  integration  with 

respect  to  the  height  from  surface  of  the  earth  to  the  height  of 

the  satellite  h^.  The  elevation  p  is  determined  (4)  .by  position  of 

the  satellite  (h^.iiQ).  The  linear  integral  I(p,T:^)  is  dependent  on 

the  coordinate  of  the  receiver  r.  and  the  elevation  of  the 

1 

satellite  p.  The  linear  Integral  is  the  complete  phase  (j>  (1); 
here,  the  reconstructed  function  P  will  be  proportional  to  N. 
Since  there  cannot  be  a  large  number  of  receivers  and  the  range  of 
angles  p  is  limited,  it  is  inadvisable  to  examine  methods  for 
analytical  inversion  of  such  linear  Integrals  and  methods  for 
integral  transfoiroations.  In  given  case  small-aspect  tomography  is 
intended  from  the  start  to  solve  the  problem  in  discrete  form  and 
to  use  algebraic  reconstruction  algorithms  or  methods  for 
expansion  into  finite  series. 

At  first  we  will  examine  the  possibility  for  replacing  ray 
(T)  by  straight  line.  The  ray  became  curved  after  the  switch  to 
the  new  coordinates  (h,^)  convenient  for  solution  of  the  of  the 
discrete  problem.  The  straight  ray  connecting  the  receiver  (h^=0, 
'Cj^=0)  and  the  satellite  (h^,  Tq)  is  defined  by  the  function  h*  (1:)= 
=  1  ctg  0Q,  which  differs  from  the  dependence 

hCT)  =  R  [  sln(j)/sln(t|>~T:/R)  -13, 

where  h(n:o)  ^  h’CT^)  =  h^.  landing  in  powers  of  the  small 
component  r/R  «  (j),  we  find  that  the  height  difference  Ah  between 
the  two  trajectories  is  expressed  by  the  formula 


^0  o 

Ah  =  h’(T:)-h(T)  «  — (“g-  +  ctg  (j)Q)T: - (-^-  +  ctg^<])Q)  +  0(-^). 

R  R  R*^ 

In  the  middle  t:=Tq/2  of  the  tra;Jectory  where  t|)Q=  'ic/4,  11^=1000  tan. 
Ah  <y60  tan.  Prom  this,  the  division  of  the  ionosphere  into  vertical 
increments  should  significantly  exceed  Ah  if  the  curvature  of  the 
ray  in  coordinates  (h,T)  is  considered  or,  otherwise,  if 
"curvature"  of  the  polar  coordinates  In  the  region  of 
reconstruction  using  a  straight  ray  is  not  considered.  The  total 
phase  will  be  even  more  sensitive  to  a  failure  to  consider  the 
curvature  of  the  ray  if  we  attempt  to  reconstruct  N  from  the 
complete  phase,  since  the  lengths  of  the  curved  and  straight  rays 
will  differ  significantly.  In  short,  consideration  of  the 
curvature  of  the  polar  coordinates  in  the  region  of  reconstruction 
or  of  the  curvature  of  the  ray  in  coordinates  (h,*!;)  in  ray 
ionospheric  RT  of  global  structures  is  necessary,  which  is, 
unfortunately,  not  taken  into  account  in  a  number  of  works 
t13,17,18]. 
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1  •  ISie  Bolution  of  the  direct  problem  of  radio  wave 
propagation  for  the  Ionospheric  Tomogjraphy. 


The  purpose  of  this  section  is  to  describe  the  technique  and 
program  for  solving  the  direct  problem  of  ionosphere  radio 
probing,  i.e.  the  problem  of  obtaining  the  phase  and  the  phase 
difference  from  the  date  on  the  electron  density  distribution.  In 
view  of  calculation  the  direct  problem  solution  amounts  to 
calculatiry  the  integral  (9),  i.e.  the  pliase  or  the  difference  of 
such  integi'als  (the  phase  difference)  for  a  small  variation  of  the 
satellite  position  angle  p.  For  the  computer  modeling  of  BT 
problems  it  is  necessai'y  to  calculate  a  series  of  such  integrals 
for  arbitrary  positions  of  the  receivers  and  the 
satellite, therefore  we  accomplished  the  progi’am  for  calculating 
the  phase  and  the  phase  difference  for  arbitrary  positions  of  the 
receivers  and  the  satellite. 

Since  the  integrand  lias  no  peculiarities,  the  calculation  can 


be  performed  using  the  rectangle  technique  or  the  Simpson  method. 
Let  us  evaluate  the  necessary  step  of  numerical  integration  and 
the  accuracies  obtained  in  this  way.  It  is  well  Known  that  errors 
of  numerical  integration  by  the  rectangle  method  and  the 
Simpson  method  are  equal  to: 


^0  P  ^ 

~  rrS.  * 

12  m 

where  h^  is  the  integration 
height),  m  is  the  number  of  the 


h, 


D  F 


(4) 


= 


(10) 


~  28S0  * 

interval  length  (the  satellite 
integral  discretisation  elements. 


Integration  errors  made  using  the  rectangle  and  Simpson  methods 
are  proportional  to  the  values  of  the  second  and  fourth 


derivatives,  respectively,  of  the  integrand  at  a  certain  point 
within  the  integration  Inter’^al,  For  an  approximate  estimation  of 


errors  It  Is  sufficient  to  evaluate  the  second  and  fourth 
derivatives  as  a  result  of  dividing  a  characteristic  value  of  the 
function  F  hy  the  square  or  the  fourth  power  "a",  respectively,  of 
the  structural  irregularities  of  F.  Limitations  associated  with 
diffraction  effects  as  well  as  those  of  the  linear  tomography 
problem  make  it  impossible  to  reconstmct  details  smaller  than 
10-20  tan  using  the  method  of  ray  RT  [1,23,  therefore  a  ^  10-20  km. 
The  value  of  m  is  equal  to  the  result  of  dividing  h^  by  the 
integration  step  Ah.  Hense  the  estimations  of  absolute  and 
relative  errors  are: 


s  iS -  (Ah/a)^, 

r  12 


h^P 


^  (Ah/a) VI 2, 


e_  s  -  (Ah/a)’. 


B 


2880 


(11) 


6  S 


v 


«  (Ah/a)’/2880, 


When  the  Integration  step  Ah=0.5  km  and  a=20km  the  relative  error 
of  the  rectangle  method  is  0.5*10“^  (for  the  Simpson  method  it  is 
smaller  than  10“^^)  which  is  quite  satisfactory  for  RT 
applications. 

For  further  computer  modeling  of  the  RT  problems  it  is 
necessary  to  use  a  set  of  appropriate  electron  density 
distributions  models  of  the  lonoshere.  Naturally,  it  is  impossible 
to  make  up  a  complete  set  of  all  the  clncelvable  ionospheric 
conditions,  and  this  study  was  not  aimed  to  do  so.  For  our  purpose 
-  to  illustrate  the  applications,  which  Includes  the  main 
structural  features  (a  trough,  localized  natural  and  man-made 
irregularities  and  groups  of  irregularities).  Horewer,  the 
presented  package  of  programs  makes  it  easy  to  design  many  other 
structural  types  and  to  extend  this  "ZOO”  as  far  as  possible  using 
the  available  "details".  In  what  follows  the  available  "details" 
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and  the  models  used  are  described. 

Description  of  the  set  of  the  models  used  for  reconstruction. 

The  "parabola"  (with  a  discontinuity  of  the  first  derivative) 
and  the  "cosine  square"  (with  no  discontinuity  of  the  first 
derivative),  a  "gaussian"  are  used  as  the  functions  to  describe 
localized  irregularities.  Gross -sect ions  of  the  constant  value  for 
such  irregularities  may  be  arbitrary  oriented  ellipses. 

No.1.  A  simple  model  of  the  ionosphere  with  a  troi:igh  and  a 
positive  irregularity  at  the  left  edge  of  the  trough  is 
represented  in  Pig.1  in  isolines  in  the  lO^cm”^  units. 

No. 2.  A  model  of  the  ionosphere  with  a  trough  and  three 
Inhomogenelties  (all  of  them  beiiig  of  the  "cosine-square"  type, 
the  positive  one  is  at  the  left  edge  of  the  trough,  the  other 
positive  one  is  also  at  the  height  of  the  maximum  at  the  opposite 
edge  of  the  trough,  the  negative  irregularitity  is  located  higher 
(h=500km,  T  =  SOOkm))  representated  in  Pig. 2  in  isolines  in  the 
10®cffi~^  units. 

No. 3.  A  model  of  a  homogeneous  smooth  ionosphere  with  a  chain 
of  irregularities  in  the  region  of  the  main  maximum.  The  model  is 
represented  in  Pig.3  in  isollnes  in  the  lO^cra”*^  units. 

No. 4.  The  model  of  the  heated  ionospheric  lens  1261.  The 
maximum  of  AN  is  equal  AN=10'^cm“^.  N(h,T:)=NQ(hQ/G(h))^(1-T^/Q^eQ), 
%  ^  G(h)eQ,  G(h)=(h^+7(h-h|))'’^^,  hQ=100  tan,  7=3,  6^=0, 15  (Pig. 4). 

No. 5.  A  model  of  a  localized  irregularity  described  by  3 
"gaussians".  The  sizes  of  the  disturbed  region  are  100  x  100  tan. 

No. 6.  A  group  of  isolated  irregularities  (Pig.6)  described  by 
the  "cosine  square"  functions. 

No. 7.  A  model  (Plg.T)  of  a  localized  simple  irregularity 
described  by  two  "cosine  square".  The  sizes  of  the  disturbed 
region  are  100  x  100  tan,  the  hlght  is  200taii. 


Description  of  the  program 

for  calculation  the  phase  and  the  phase-difference  for 
arbitrary  positions  of  the  receivers  and  the  satellite. 


System  Requirements 


-  Computer: 

-  Operating  System: 

-  Memory: 


-  Hard  Disk  Space: 

-  Software: 


IBM  AT  or  compatible  (with  coprocessor) 
MS-DOS  or  PG-DOS  version  3.0  and  later 
at  least  Extended  memory  16  Mbytes 
(depends  on  geometry  and  type  of 
approximation  reconstructed  function) 
16  Mbytes 

compiler  1 .4e  and  linker  2. 2d 
NDP-P0RTRAN-386(c)  MicroWay  or  later 


1.  Program  <lntegr.for> 

This  program  solves  direct  problem,  namely,  determines  the 
model  structure  and  calculates  the  phases  or  TEC  for  phase  RT  or 
doppler  for  phase- difference  RT  on  the  model  structure. 

Input  parameters  and  files: 

Answer=1  -  Phase  RT:  Tec_ph=1  ~>  calculation  of  TEC 

Tec_ph=2  — >  calculation  of  Phase 
Answer=2  -  Phase-difference  RT 
Answer,  Tecjph  are  introduced  from  the  screen 
file  <name_P.int>  -  names  of  output  files  for  integrals  for  each 

receiver 

NREC  -  number  of  receivers 
KRAY  -  max  number  of  rays 
RAY  -  array  of  rays  for  each  receiver 
NP  -  number  of  discrets  on  the  horizontal  grid 
NR  -  number  of  discrets  on  the  vertical  grid 
NirMAX  -  max  number  of  local  irregularities  for  model 
PP  -  array  of  receivers ’s  polar  angles  in  degrees 
(  PP(1)=0.  ->  for  the  first  receiver) 

RZ  -  Earth's  radius  in  km 


RO  -  satellite  orbit  in  polar  coordinate  system  in  km 
HPISTjHPIN  -  Initial  and  final  altitude  of  reconstructed 
structure  in  km 

P1  -  horizontal  distance  to  the  left  from  first  recelYer 
in  km  (P1<0) 

P2  -  horizontal  distance  to  the  right  from  last  recelYer 
In  km  (P2>0) 

Dconst  -  the  Yalue  for  deteimination  of  doppler  (shift  of  rays) 
NJ  -  number  of  discrets  on  the  one  ray  for  calculation 
of  Integrals 

RMAX  -  altitude  of  max  electron  density  in  km 
The  model  structure  is  determined  by  function  <FMODEL>  which 
uses  functions  <PUNG1>,  <PMC2>,  <HOMPAR>  or  <HOMGOS>  and  input 
file  <P_mod>. 

Input  file  <F_mod>  contains  parameters  of  model,  for  example: 

«  ’number  of  irregularities  and  Zmax  -  max  Yalue  of  parabol’ 

1  1. 

'local  irregularities:  Hirreg,  Tirreg,  A,  B,  Z' 

360.  36.  100.  100.  0.08 
'trough  — >  TGFJ  ,  TGE2  in  km' 

-20.  640.  » 
line  2  in  file:  1  1 . 

Nirreg  (number  of  Irregularities )=1 
Zmax=1 . 

line  4  in  file:  350.  36.  100.  100.  0.08  - 
Hir  (altitude  of  irregularity  in  km) =350. 

Tir  (horizontal  coordinate  of  irregularity  in  km)  =36. 

Alr=100.-  Yertical  size  of  irregularity 
Bir=100.-  horizontal  size  of  irregularity 
Zmax=0.08  -  max  Yalue  of  Irregularity 
line  6  in  file 
-20.  640. 

TGR1=-20.,  TGR2=640  -  horizontal  coordinates  of  trough  in  km 

RMAX,  RM,  B1,  B2,  ZSM  -  parameters  for  function  <FUN01> 

(model  on  h) 

DEPTH,  G0NS1 ,  C01JS2  -  parameters  for  function  <FUNG2>  (trough) 
DEPTH  -  depth  of  trough 
G0NS1,00NS2-  the  egdes  of  trough 


<HOMPAR>  -  approximation  of  local  irregularity  by  parabol 
<HOMCOS>  -  approximation  of  local  irregularity  by  cos^ 

The  combinations  of  <PM01>,<FUNG2>,<HOMPAH>,<HOMCOS>  make  it 
possible  to  obtaine  different  model  structures. 

Subroutine  <DEFPSI>  determinates  the  angles  of  rays  on  the 
satellite’s  coordinates 

Subroutine  <DEPINT>  calculates  either  Phase  or  TEG  or  Doppler  and 
uses  subroutine  <INTERG>. 

Subroutine  <INTERG>  calculates  integral  only  one  ray  on  the  model. 

Output  files: 

MOD.GRD  -  file  of  model  structure 

Files  <Fint>  -  arrays  of  either  Phases  or  TEG  or  Doppler 

Execution 

f77  Integr.for 

RUM 

Integr.exp 


2.  DesigQ  of  the  different  versions  of 
the  KP  operators  (nia trices). 

We  shall  consider  the  problems  of  the  design  of  the  different 
versions  of  the  RT  operators  (matrices),  beginning  with 
discretization  procedure  for  equations  (1).  We  perform 
digitization  of  the  linear  Integrals  I(p.‘ij)  (9)  according  to  the 
position  of  the  satellite,  which  Is  dependent  on  the  coordinates 
angle  'tQ.?/R.  The  set  of  coordinates  of  the 

satellite  Tqj  Is  recalculated  according  to  (4)  Into  a  series  of 
elevation  p  .  .  of  the  satellite  from  receiver  1: 

X  J 

(R  +  Iiq)  cos(a^  -  -  R 

tg  p  =  -  ,  (12) 

(R  +hQ)  sln(a^  -  a^^.) 

The  sets  of  elevation  of  all  receivers  define  a  series  of  discrete 
values  of  the  linear  Integrals  s  ’i^).  The  simplest 

method  for  digitization  of  the  sought  function  In  a  fixed 

rectangular  (m^  x  ng)  grid  is  to  replace  it  by  a 
plecewlse-oonstant  approximation,  or  to  represent  P  by  a  system  of 
(mg  X  Hq)  basis  functions  equal  to  unity  in  certain  rectangle  and 
zero  in  all  others.  The  rectangular  reconstruction  region  is 
divided  Into  heights  (m  ^  m^)  and  Hq  horizontal  samples  (1  <  n 
^  Hg).  Let  the  value  of  the  function  P(h,T:)  in  a  fixed  (m  x  n) 
rectangle  be  The  point  In  the  rectangles  at  which  the  samples 
of  P(h,T)  are  selected  Is  not . especially  Important;  this  may  be  at 
the  middle  of  the  rectangles  or  nodes  of  the  grid. 

The  problem  of  tomographic  reconstruction  according  to  linear 
Integrals  Is  to  determine  the  set  of  discrete  samples  CP^^}  In  the 
known  grid  according  to  the  set  Desiguatlng  the  length  of 

ray  (1,J)  In  cell  (m,n)  as  we  obtain  the  system  of  linear 

X,  J 


if 


equations 

m,n 

L  „  =  I.  ,• 

•  j  m^n  ijD 
•*-»*} 


or 


M 

* 


(13) 


Here,  "renumbering”  of  the  ray  (i,3)-*  J  and  the  sells  of  the 
ionosphere  (m,n)-^  M  is  performed  in  the  second  equation.  The 
repeating  indices  are  understood  as  summation.  The  number  of  rays 
is  determined  by  the  parameters  of  the  recording  system.  The 
coefficients  are  calculated  according  to  the  given  rays  and 
cells  into  which  the  ionosphere  is  divided.  System  (13)  may  be 
either  overdeteimined  or  sub-definite. 

Thus  if  there  is  a  possibility  to  determine  the  linear 
integrals  (1),  the  problems  of  the  ray  RT  are  reduced  with  the 
help  of  discretization  procedure  to  solving  systems  of  linear 
equations.  But  the  problem  of  ionospheric  RT  according  to 
phase-difference  or  Doppler  measurements  cannot  be  solved  by  such 
scheme  with  a  piecewlse-constant  approximation.  The  fact  is  that 
the  data  here  will  be  derivatives  of  linear  integrals  of  type  (9): 
D  =  dl/dUQ,  or  finite-difference  ratios  of  the  increment  AI  of  the 
linear  integrals  to  the  increment  Aa^  of  the  satellite  coordinate. 
The  Doppler  frequency  =  d4>/dt  measured  in  the  experiment  is 
determined  by  the  phase  derivative  (1).  The  relation  between  the 
angle  of  a  satellite  moving  uniformly  along  a  circular  orbit 
with  velocity  Vq  and  the  time  Uq  =  VQt/(R+hQ)  makes  it  possible  to 
express  the  Doppler  frequency  Q  by  means  of  the  derivative  with 
respect  to  the  angle  of  the  satellite 

^0  d(}) 

R+h^  ddg’ 

from  which  these  phase -difference  tomography  data  are  proportional 
to  AI/AUq.  The  derivatives  of  the  linear  Integrals  in  a 


plecewlse-constant  approximation  of  the  sought  function  F  will  he 
discontinuous.  This  Is  because  each  linear  Integral  is  the  sum  of 
integrals  OYer  the  set  of  cells.  As  the  elevation  of  the  satellite 
changes,  the  ray  encounters  a  new  cell;  the  integral  with  respect 
to  unity  of  this  cell  is  a  continuous  function  of  the  angle  of  the 
satellite  a^,  but  the  derivative  of  the  linear  integral  with 
respect  to  will  contain  a  discontinuity  when  the  ray  contacts 
the  comer  of  each  cell.  Therefore,  the  plecewlse-constant 
representation  of  the  function  to  be  reconstraoted  does  not  make 
it  possible  to  analyse  the  phase-difference  problem. 

Phase-difference  measurements  require  ■  higher-order 
interpolation  than  the  plecewlse-constant  representation  of  the 
function  to  be  recorded-.  Correspondingly,  the  matrix  Ljjj  for 
transition  from  the  function  to  be  reconstructed  to  linear 
Integrals  should  be  calcsilated  differently  so  as  to  ensure 
continuity  of  linear  Integrals  with  respect  to  the  coordinate  of 
the  satellite  (or  elevation  p).  If  the  matrix  of  the  direct 
problem  Ij^:  Ij  is  continuous  with  respect  to  the  angle  of 
the  satellite  then  in  place  system  (13)  it  is  possible  to 
obtain  a  system  for  phase-difference  or  Doppler  data  by 
differentiating  (13)  with  respect  to  the  angle  Uq: 

Here,  D^  s  AIj/Aa^  are  Doppler  data  and  s  Al^/AaQ  is  finite 
-difference  ratio  (or  derivative)  of  the  matrix  to  the 
increment  of  the  angle.  The  Doppler  data  are  determined  not  only 
by  the  change  in  the  complete  phase  related  to  integral  electron 
concentration  along  the  ray,  but  also  the  local  electron 
concentration  at  the  point  of  the  satellite.  The  correction  to 
the  Doppler  data  is  equal  to  the  product  of  times  the  velocity 


component  of  the  satellite  directed  along  the  ray  -  X  t  H  cos  (a.  + 
+  p  -  a^).  This  correction  can  be  Inserted  into  the  iteration 
algorithm  and  the  values  of  at  the  boundary  h=hQ  of  the 
ionosphere  obtained  in  the  iteration  process  will  constantly 
"correct"  the  measured  Doppler  frequency  values  [253 .  The  question 
of  whether  such  a  correction  of  the  data  should  be  made  can  be 
answered  depending  on  whether  we  use  direct  measurerQents  of  the 
Doppler  frequency  or  the  phase  difference  derived  from  ptiase 
measurements.  In  the  former  case  such  a  correction  should  be  made, 
in  the  latter  -  not.  Here  such  a  correction  shall  not  be 
performed;  horewer,  it  would  not  be  difficult  to  Introduce  it  when 
using  direct  doppler  meas’orements , 

In  what  follows  we  shall  briefly  describe  the  methods  of 
constructing  the  operators  that  are  smooth  by  the  satellite  angle; 
here  are  the  examples  of  cons  tract  ing  the  matrices  of  the 
transition  from  the  recvonstruoted  function  to  the  linear  integrals 
(matrices  of  projection  operators)  useful  as  for  the 
phase-difference  RT  and  also  for  the  phase  RT.  In  this  section  we 
consider  examples  of  constructing  smooth  projection  operators  of 
the  direct  problem.  One  must  construct  such  a  L-r,,  matrix  that 
would  provide  the  smooth  of  lliaear  integral  over  the  satellite 
angle.  In  the  beginning  section  we  consider  contribution  on  the 
basis  of  triangular  elements,  at  the  end  of  the  section  other 
possible  variants  be  outlined.  We  will  proceed  to  calculation  of 
the  matrix  of  the  difference  problem,  which,  as  was  already 
noted,  should  be  detenulned  from  the  increment  of  the  matrix  L 

JM. 

Which  is  smooth  with  respect  to  the  angle  of  the  satellite. 
Smoothness  of  the  matrix  can  be  ensured  by  introducing  finite 
triangular  elements  for  representation  of  the  function  F(h,T:), 
i.e.,  when  the  sou^t  function  is  replaced  by  a  plecewise-planar 


approximation.  The  smooth  function  P(h,T)  is  replaced  by  a 
contlnuouB  polyhedTal  approximation  surface,  according  to  which 
the  derlvatlYe  with  respect  to  the  satellite  angle  of  the  linear 
Integi'^als  is  already  a  continuous  fimction.  Triangular  elements 
are  obtained  naturally  from  a  grid  of  rectangles  by  dlYiding  each 
of  them  in  half  diagonally.  The  function  P(h,T)  within  each 
triangular  element  la  replaced  by  linear  approximation 

P{h,T:)  =  a  +  bn:  +  ch.  (15) 


The  Yalues  of  the  coefficients  (a,  b,  c)  in  each  finite  element 
are  determined  from  system  of  three  equations  for  three  boundary 
points.  It  is  simple  to  immediately  write  expressions  for  the 
coefficients  in  a  glYen  finite  element.  These  expressions  differ 
slightly  for  triangilar  elements  of  two  types;  those  occupying 
cells  "below”  or  "above".  We  stipulate  that  cell  (m,n) 

is  divided  by  a  diagonal  )  •-  ,  h^^)  running  downward 
and  left  -to-right,  into  two  triangular  elements:  the  "lower"  and 
"upper"  elements.  Then,  in  the  lower  (m,n)  element 


P(li,i:)=  P 


m,n 


P  -  -  P  P  -  P 

ra+1  ,n  m,n  ra,n+1  m,n 

+ - - (T-m  )  +  — 


m' 


Ah 


r-(h-h^)  (16) 


and  in  the  upper  (m,n)  element 
P(h,T;)  =  P 


p  —  p 

ra+1  ,n+1  in,n+1 


ra+l  ,n+1 


(''-Vi )  + 


+ 


F  —  V 

m+1  ,n+1  m+1 ,11 

Ah  ~ 


^^“^n+1 ^  * 


(IT) 


As  before,  to  simplify  the  notation  we  will  renumber  the  values  of 

the  samples  below:  P„  p„, 

ra,n  M’ 


(ffi+1,n)  -*  (M+l),  (m,n+1)  (M+AM), 


i9 


(ffi+1,n+1)  (M+AM+1),  where  AM  is  the  number  of  cells  horizontally 

In  one  row. 

The  linear  Integral  Ij  (9)  Is  the  sum  of  the  Integrals  with 
respect  to  all  finite  elements  which  Intersect  ray  J: 


Ij  =  f  7(h)  PCh.T)  dh, 

where  'r(h)=(R+h)tR2sln2p+2Eh+h2r''^?  while  F(h,t>  Is  represented 
m  the  form  of  plecewlse-planar  approximations  (16), (17)  In  each 
finite  element.  The  result  of  Integration  of  such  an  approximation 
In  lower  element  M 

J  7(h)P  dh  =  Jq  Fa  +  Jx(%i  - 
and.  In  upper  triangular  element  M 

J  7(h)I'  dh  =  ^M+AM+1  ^^M+AM+1  "  ^M+AM^ 


*^h  ^^M+AM+1  ^M+1^ 

Here,  J^,  J^,  following  Integrals: 


h 


\l+1 


J.  =  I-  J  TW  ['i:(h)-\,ldh;  ^  J  7(h)['t(h)-T^l]dh; 

'  £x  V 


At: 


h 


h 


1  -  1 
=  -  J  7(h)  (h-hj^idh:  -^h  =  -f 


(20) 


Integration  with  respect  to  a  lower  finite  element  begins  at  the 
lower  boundary  of  the  cell  h=hj^  and.  ends  at  height  h,  where  the 
ray  leaves  the  lower  element.  Integration  with  respect  to  the 
upper  finite  elements  begins  at  this  height  and  ends  at  the  height 
of  the  upper  boundary  of  the  cell.  We  will  recall  that  7(h)  and 
all  Integrals  with  respect  to  cell  M  are  functions  of  the 


10 


elevation  p,  or  the  number  of  the  ray  J. 

After  Integration  (18)  with  respect  to  ray  J  In  the  lower 
element  M  the  value  (Jq  -  -  Jj^)  is  entered  into  the  coefficient 
Ljjj,  since  it  Is  a  coefficient  for  Fjj.  Correspondingly,  is 
entered  into  Lj  and  Jj^  into  Lj  However,  in  integration 
with  respect  to  one  lower  finite  element  M  these  coefficients  L 
are  still  not  completely  determined.  It  is  easy  to  understand  that 
each  sample  falls  into  three  lower  and  three  upper  adjacent 
finite  elements.  Only  after  integration  with  respect  to  ray  J  in 
all  finite  elements  is  it  possible  to  conqpletely  form  the 
coefficient  Ljjj  from  six  around  where  the  ray  fell.  Integration 
with  respect  to  the  upper  finite  element  M  with  respect  to  ray  J 
(19)  makes  the  contribution  (J^  +  +  J^)  to  the  coefficient 
Lj  jf+jjju+i ;  the  contribution  (-J,J-)  to  the  coefficient  Ij  and 
the  contribution  (-J^)  to  the  coefficient  Lj  .  The  Integrals 
with  respect  to  all  rays  of  type  (20)  can  be  calculated  by  various 
numerical  methods;  in  view  of  the  smoothness  of  7(h)  and  the 
piecewlse-planar  approximation  of  F,  it  is  sufficient  to  use  the 
trapezoid  or  Simpson  method.  Here,  in  each  integration  step  Ah  it 
is  necessary  to  verify  that  the  ray  does  not  exceed  the  limits  of 
the  finite  element. 

Performing  numerical  integration  with  respect  to  all  rays,  we 
obtain  the  matrix  Ljjj.  The  matrix  is  related  to  the  set  Ca^) 
of  positions  of  the  satellite  and  the  corresponding  series  of 
rays.  It  is  also  possible  to  calculate  the  matrix  I'  for  another 
set  of  close  positions  of  the  satellite  with  a  fixed  Increment 
{aQ+Aag}.  After  this,  we  determine  the  matrix  for  phase-difference 
tomography  problem  Ajjj  =  (L^  -  Ljjj)/AaQ. 

The  projection  operator  or  Ljj^  matrix  can  also  be  built  on 
the  basis  of  approximations  of  higher  order  than  that  of  (15).  For 


example,  one  can  use  a  two-dimensional  approximation  in  the  form 
of  the  product  of  linear  functions  or  the  product  of  cubic 
splines.  Then  the  function  I'(h,'i:)  takes  the  form 

m,n=0 

Inside  an  arbitrary  (m,n)  rectangular,  representing  the  function 
through  the  normalized  coordinates  x=(t:-t:jjj)/At:,  y=(h-hj^)/Ah,  we 
can  obtain  the  following  representation  through  the  values  of  the 
function  at  the  four  angular  points  {x,y)=  {(0,0);  (0,1);  (1,0); 


(1,1)}. 


F(x,y)  =  Poo^oo+ 


^00^00'^ 


y  y 
^oo^oo'*’ 


xy  zy 
^00^00'^  ^ 


Here  the  subscrips  refer  to  the  coordinates  (x,y)  of  the  angular 
points,  F  is  the  value  of  the  function  in  corresponding  point,  F^, 
F^  are  the  values  of  the  function  partial  derivatives  with  respect 
to  x,y,  F^  is  the  value  of  the  function  partial  derivative  of  the 
second  order  with  respect  to  x,y.  The  total  sum  (21)  will  contain 
16  summands,  P3Lp(3:,y)  are  corresponding  polynomials  of  the  power 
not  exceeding  3.  We  shall  not  write  the  mentioned  polynlmlals 
completely,  four  examples  will  be  sirffice  for  illustration: 


Pqq=  6x^y^-  6x^y^+  2y^-  3y^+2x^-3x^+  1 , 

Pqq=  ExV-  4xy-  3x^y^+  6x^y^+  2xy^-  3xy^+x^-2x^+  x, 

Pqq=  2xV-  Sx^y^-  4xy+  6x2y2+  2x^y  _  +y3_2y2^.  y^ 

p^=  x^y^-  2x^y^-  2x^y^+  4x^y^+  x^y  -  2x^y  +  xy^-  2x  y^  +  xy. 

Then,  integrating  in  each  cell  of  the  given  polynomials  we 
can  p2X)duce  the  corresponding  elements  of  the  matrix,  as  In  the 
case  (16-20).  Note,  that  now  it  is  not  only  the  values  of  the 
function  F,  but  the  values  of  the  mensloned  derivatives  that  are 
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imknown,  i.e.  such  representation  makes  it  possible  to  find  the 
function  and  its  first  derivatives.  The  matrix  for  the  product  of 
linear  approximations  can  be  constructed  In  a  similar,  even 
simpler,  way,  this  being  In  fact  a  particular  case  of  that 
described  above,  where  the  summation  In  the  formula  la  made  up  to 
1  rather  than  up  to  3. 

Thus,  the  following  operators  (matrices)  for  solving  the  RT 
problems  are  described: 

A  -  Is  built  with  the  piece-constant  approximation, 

B  -  Is  built  with  the  piece-planar  approximation, 

-  Is  built  with  the  linear  product  approximation, 

D  -  Is  built  with  the  cubic  spline  product  approximation. 

Pro;Jectlon  operators  with  approximations  of  higher  orders 
allow  a  better  approximation  of  the  operator  of  the  direct 
problem,  i.e.  they  enable  us  to  come  closer  to  the  true  operator 
of  the  direct  problem.  Tables  1,2  shows  examples  of  calculations 
of  the  direct  problem  for  the  model  2  and  5  with  the  help  of 
different  operators:  A,  B,  G,  D.  Errors  of  the  numerical 
simulation  can  appropriately  be  characterized  by  numbei^  6,  which 
shows  the  deviation  of  the  function  being  reconstructed  E  from  the 
true  function  P:  6=jE-f|/fE|.  The  norms  of  the  spaces  1^  and  l“ 
can  be  helpfully  used  (dg  s  d(l^)  and  s  6(l“)).  We  use  the  data 
(hQ=1 000km)  from  three  receivers  with  coordinates  0km,  ^2=  4T5 
km,  1435  km  (for  model  2,  this  geometry  is  similar  to  the 
geometry  of  Munnansk-Moscow  RT  experiments  t1,2])  and  1:^=0, 
't2=240,  0:2=480  (model  5).  Prom  the  results  given  In  this  tables 
one  can  clearly  see  the  increased  accuracy  In  solving  the  direct 
problem  for  operators  with  higher  oders  of  approximation. 

One  can  see  that  the  transition  to  higher  orders  makes  it 
possible  to  significantly  Improve  the  solution  of  the  direct 
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problem.  However,  as  the  approximation  order  Increases  the  matrix 
becomes  more  complicated  and  less  rarefied,  which  can  Impair  the 
solution  of  the  inverse  problem.  One  cannot  say  in  advance  which 
operator  will  be  more  preferable  when  solving  the  inverse  problem, 
because  at  the  beginning  of  the  Increase  of  the  approximation 
order  the  fmction  is  approximated  better  but  the  matrix 
properties  for  solving  the  Inverse  problem  become  worse.  Operators 
must  be  chosen  by  means  of  computer  simulation  to  be  Illustrated 
in  the  next  section. 


Description  of  the  program  for  calculation 
of  the  different  versions  of  the  RT  matrices  (operators). 

I  2.  Program  <roatric .  f or> 

/  This  program  desigis  different  versions  of  the  operators 
j  (matrices)  for  phase  RT  or  phase-difference  RT. 

I 

1 

'  Input  paiumeters  and  files  s 

Parameters  <  Answer,  JfF,  NR,  NREO,  NJ,  RAY,  HPIST,  HPIN,  P1 ,  P2, 
R2,  RO,  Dconst>  are  similar  to  same  parameters  of  program 
<INTEGR.POR>. 

Pile  <name_F.mat>  contains  names  of  output  files 
MOD.GRD  -  input  file  of  model  structure  (program  <INTEGR.POR>) 
APTYPE  -  type  of  approximation  of  reconstructed  structiire,  it  is 
introduced  from  the  screen 

Subroutine  <DEPPSI>  determinates  the  angles  of  rays  on  the 
satellite's  coordinates. 

Subroutine  <DEPMAT>  determinates  the 

matric  with  corresponding  approximation  for  one  receiver 


Output  files I 

Files  <F_matr>  -  arrays  of  matric  of  corresponding  approximation 

for  each  receiver 

File  <Fparam>  -  parameters  of  matrices  for  each  receiver 
Files  <F_lnt>  -  results  of  multiplications  of  calculated  matric 
and  model  structure 

File  <F_st>  -  array  (LST)  of  number  of  all  rays  which  cross  the 
corresponding  dlscret  of  model  (SIRT  algorithm) 

Execution 

f77  matric. for  raylnt.ob;! 

RUN 

matric. exp 


3.  The  solution  of  the  inverse  problems  of  the  phase  and 
phase-difference  RT. 

It  should  be  noted  that  the  earlier  authors  solved  problems 
of  tomography  using  linear  integrals,  which  amounted  to  solving 
systems  of  linear  equations.  In  this  respect,  therefore,  the 
problem  of  pure  phase  RT  is  entirely  equivalent  to  the  known 
problems  on  linear  integrals.  The  peculiarity  of  ionospheric 
upplications  is  revealed  in  the  nature  of  possible  errors  by  a 
conctant  in  determining  the  phase,  which  will  be  dealt  with  in  the 
next  section. 

In  this  section  varios  methods  of  solving  Inverse  problems  of 
RT  will  be  brlfly  described  and  analysed.  The  results  obtained 
using  various  reconstruction  methods  and  the  sensitivities  of 
these  methods  will  be  compared,  which  is  of  Interest  for 
ionospheric  application.  Solution  of  (13-14)  in  ray  RT  of  the 
ionosphere  is  difficult  in  computational  respects.  When 
reconstructing  global  structures  with  dimensions  of  thousands 
kilometers  and  a  sampling  interval  of  tens  kilometers,  the 
matrices  of  such  systems  contain  up  to  10^-10*^  elements,  but  are 
rather  empty.  There  are  a  significant  number  of  both  direct  and 
Iterative  methods  for  solving  system  of  linear  equations  like 
(13-14).  Nevertheless,  Intensive  development  of  the  theory  and 
practical  methods  and  algorithms  for  solution  of  these  systems 
continues  at  present.  There  are  a  large  number  of  diverse 
iterative  methods  for  solving  systems  of  linear  equations.  Many  of 
them  have  been  tested  in  ray  RT  problems.  As  was  previously  noted, 
tomographic  methods  have  been  developed  most  intensively  in 
seismlcs,  where  various  methods  for  solution  of  these  system  have 
been  used.  Here,  we  can  cite  algebraic  reconstruction  technique 
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(ART),  including  with  relaxation  and  for  systems  of  inequalities; 
siimiltaneous  iteration  reconstruction  technique  (SIRT); 
multiplicative  algebraic  reconstruction  technique  (MART); 
block-iterative  algorithms,  reconstruction  on  the  basis  of  the 
Bayes  approach,  algorithms  for  regularization  of  the  mean  square 
error,  algorithms  for  optimization  of  image  entropy,  etc.  E363. 
Good  practical  results  have  been  obtained  in  seismic  tomography 
using  methods  for  minimization  of  the  iterative  corrections  in 
various  matrices  with  variants  of  weighted  correction,  ray 
weigh-ing  and  inter-iteration  smoothing  13].  Investigations  have 
been  made  of  the  characteristics  of  the  spectra  of  the  matrices, 
the  resolution  of  the  reconstruction  systems  and  the  miqueness  of 
the  reconstruction  possible  ”at  the  limit”,  i.e.,  with  an 
infinite  increase  In  the  number  of  measurements  [53. 

Since  there  are  numerous  methods  of  solving  systems  of  linear 
equations,  it  appears  to  be  impossible  to  apply  all  the  Imown 
approaches  within  one  study.  It  should  be  pointed  out  tliat  the 
main  results  of  the  investigations  perfomed  are  weakly  dependent 
on  the  applied  methods  for  solving  systems  of  linear  equations, 
here,  therefore,  we  sliall  apply  the  most  widely  known  methods: 
ART,  SIRT,  MART. 

Before  giving  examples  of  RT-reconstructlons,  it  is  to  be 
noted  that  in  some  cases  there  is  a  possibility  of  an  accurate 
determination  of  the  absolute  pliase.  It  is  possible  to  determine 
the  absolute  phase  at  Inhomogeneities  when  reconstructing  the 
structure  of  sufficiently  large  localized  artificial  (releases, 
heating,  etc.)  or  natural  Inhomogeneities  of  the  ionosphere 
arising  during  the  time  between  flights  of  the  satellite.  Such  a 
formation  being  localized  in  space  provide  the  possibility  of 
solving  the  problem  witliout  additional  a  priori  assumptions 


conceralng  the  Inhomogeneity.  By  recording  the  data  before  and 
after  the  emergence  of  the  disturbance  and  making  subtraction  of 
data,  one  can  find  the  contribution  of  the  localized  inhomogeneity 
being  reconstructed,  if  the  ionosphere  changes  little  between  the 
flights. 

Fig. 8  shows  the  reconstruction  of  a  heated  ionospheric  lens 
by  phase  RT  after  20  ART  iterations  using  the  data  from  three 
receivers  with  coordinates  t:^=0  km,  1:2=45  km,  ^^=90  km,  hg  = 
1000km.  The  height  of  the  disturbed  region  varied  from  80  to  220 
km.  The  zeroth  initial  guess  was  given.  The  size  of  the  division 
discrete  was  lOkmxIOkm.  The  reconstruction  differs  little  from  the 
model  structure,  the  relative  reconstruction  errors  being  Og  = 
=0.07,  5^  =0.09.  Also  the  simulation  of  the  reconstructions  of 
other  localized  formulations  was  carried  out,  thus  showing  the 
possibility  to  reconstruct  structures  of  localized  objects 
emerging  between  flights  of  the  satellite  using  data  on  the 
absolute  phase  of  this  object. 

In  what  follows  the  re3u3.ts  of  modeling  the  RT 
reconsotraction  using  different  operators  A,B,C,D  aru  presented. 
Tables  2  and  3  show  examples  of  calculations  of  the  reconstruction 
results  for  the  models  2,  5  and  7  with  the  help  of  different 
operators:  A,  B,  0,  D.  Errors  of  the  numerical  simulation  be 
characterized  by  number  6^  and  6^  ,  which  shows  the  ratio  of  the 
deviation  of  the  function  being  reconstructed  F  from  the  true 
function  F.  The  reconstruction  results  must  be  compared  in  the 
norms  approaching  to  the  integral  ones,  i.e.  the  norms  must  be 
used  instead  of  1^  and  instead  of  1“  (the  corresponding  numbers 
Ag  and  A^  instead  of  numbers  6g  and  In  other  words  the 
functions  must  be  compared  using  a  much  finer  grid  than  that  used 
for  reconstruction.  Otherwise,  if  the  results  are  compared  in 
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reconstruction  points  only,  we  may  come  to  erroneous  conclusions. 
We  use  the  data  from  three  receivers  (similar  to  Mnrmansk-Moscow 
KP  experiments  E1,2])  with  coordinates  0:^=  0  km,  475  Inn,  1^= 
1435  km  (for  model  2),  t:^=0,  t:2=240,  t:2=480  (model  5)  and  t^=0, 
t:2=250,  ^2=500,  1^=750,  Tg=1000,  (model  7).  Pigs. 9  (A,B,C,D)  show 
the  RT  reconstruction  results  for  different  approximations  (A,  B, 
0,  D,  respectively)  of  model  7.  Pigs. 10  (B,C)  show  the  RT 
reconstruction  results  for  different  approximations  (B,  0)  of 
model  2.  Pig.  11  show  the  RT  reconstruction  result  for 
approximation  C  of  model  5  and  fig. 12  for  approximation  (D).  Prom 
the  reconstruction  results  and  tables  2,3  one  can  see  the  increase 
in  the  accurasy  of  solving  RT  reconstruction  problems  for 
operators  with  higher  orders  of  approximation. 

The  results  of  RT  reconstruction  also  depend  on  the  applied 
approach  of  the  phase  or  phase-difference  RT.  If  the  conventional 
phase  RT  with  the  piece-constant  (A)  approximation  is  used,  the 
phase-difference  RT  with  matrices  of  the  (E,G,D)  type  has,  as  a 
rule,  an  advantage  connected  with  a  more  accurate  representation 
of  the  direct  problem  operator.  Pigs. 13  and  14  show  the  results  of 
the  RT  reconstruction  of  model  2  using  the  methods  of  the  phase 
(fig.  13)  and  phase-difference  (fig.  14)  RT.  The  homogenious 
ionosphere  model  (flg.15)  was  used  as  an  Initial  approximation. 
Phase-difference  RT  allows  a  more  precise  isolation  of  local 
extrema  and  has  a  lower  noise  level.  Note  that,  in  general,  the 
reconstruction  results  of  rather  large  structures  obtained  using 
these  methods  providing  there  were  no  errors  in  determination  of 
the  Initial  phase  are  comparable.  The  presense  of  such  an  error, 
however,  makes  the  phase  RT  method  practically  msultable,  which 
will  be  shown  by  the  results  to  be  presented  in  the  next  section. 

Let  us  analyse  the  results  obtained  by  applying  various 


methods  of  solving  systems  of  linear  eq-uations  (SLE)  for  RT 
problems.  It  has  already  been  mentioned  that  we  shall  restrict  our 
consideration  to  the  most  widely  used  techniques,  namely,  the  ART, 
SIRT  and  MART  alglrlthms.  Comparison  of  different  algorithms  of 
solving  SIiE  must  clearly  be  made  for  the  same  RT  method.  The 
results  of  comparing  the  SIS  solution  algorithms  for  the  phase  RT 
method  will  be  Illustrated  by  the  reconstruction  of  model  6.  Since 
reconstruction  errors  are  determined  both  by  the  number  of 
iterations  and  measurement  errors,  each  algorithm  should  be 
characterized  by  a  two-dimensional  error  "portrait"  providing  the 
dependence  of  a  relative  reconstruction  error  on  the  number  of 
iterations  and  the  relative  error  in  the  rlglit-hsnd  part  of 
equation.  Such  an  error  portrait  characterizes  the  convergence 
rate,  the  feasible  minimum  of  the  reconstruction  error  for 
different  levels  of  data  errors.  The  error  portraits  for  the  ART 
and  MART  algoi’ithms  are  represented  in  figs.  16,17.  One  can  see 
that  for  non-zero  errors  the  algorithms  begin  diverging  quite 
rapidly,  the  Iteration  process  must,  therefore,  be  stopped  In  the 
region  of  the  reconstruction  error  minimum,  whose  position  being 
deteiT/ilned  by  data  errors.  The  SIRT  algorithm  (Rig.  18)  is 
significantly  less  sensitive  to  data  errors  owing  to  intermediate 
averaging  of  the  results  during  the  iteration  process.  The  error 
level  in  the  SIRT,  however,  is  much  higher  {M3%)  than  those  of 
the  ART  and  MART  algorithms.  The  numerical  experiments  performed 
with  other  models  also  showed  that  the  SIRT  algorithm  is 
practically  unsuitable  due  to  the  high  level  of  reconstruction 
errors.  This  algorithm  would  be  appropriate  for  use  in  the  case  of 
severe  data  errors  (^-10%).  In  this  case,  however,  the  SIRT 
algorithm  would  reconstruct  a  highly  "averaged  obiJect"  which  bears 
little  resemblance  to  the  original  structure.  In  general,  the 
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reconstruction  results  obtained  using  the  ART  and  MART  algorithms 
are  comparable:  in  some  cases  the  ART  has  certain  advantages,  in 
other  cases  the  MART  is  more  appropriate. 

Another  advantage  of  phase-difference  tomography  over  phase 
tomography  should  be  noted,  that  is  a  higher  sensitivity  of  the 
former.  Doppler  data  are  more  sensitive  to  small  inhomogeneltles 
which  have  little  affect  on  the  phase.  Por  example,  when  the  ray 
scans  the  Inhomogenelty  AN  of  the  size  a,  the  phase  changes  by  A® 
Ar^ANa,  here  the  relative  change  of  the  phase  A®/®  ANa/Ni.'^ 

u 

"'ANa/N^L,  where  is  the  TEO  along  the  ray,  is  the  value  at 
the  maximum  electron  concentration,  L  is  the  ray  lenght 
characteristic.  As  a  result,  relative  variations  of  the  phase  are 
proportional  to  the  ratio  of  the  TEG  of  the  inhomogenelty  to  the 
TEC  of  the  whole  ionosphere.  One  should  not  expect  the  methods  of 
solution  of  (13,14)  to  be  more  sensitive  to  clianges  of  the 
right-hand  part  than  a  few  percent.  Therefore  phase  methods  would 
not  distinguish  even  sufficiently  strong  AN/Nj^"  0.1 
inhomogeneltles  with  the  size  a<L/10  (<si00kffl),  since  they  produce 
only  1%  of  phase  variations.  It  is  not  accidental,  in  our  opinion, 
that  In  the  reported  reconstructions  using  the  phase  methods 
[18-203  details  with  dimensions  of  less  than  a  few  hungred 
kilometers  are  not  revealed.  This  is  not  the  case  for 
phase-difference  measurements.  Here  total  Doppler  variations  are 
proportional  to  dc|)/dt  Ar^NL/ (L/u^ ) ,  and  Doppler  variations  at  an 
Inhomogenelty  are  "  Ar^ANa/Ca/Ug).  Then  relative  Doppler 
variations  (phase  differences)  are  proportional  to  the  ratio  of 
electron  concentrations  Thus,  the  phase-difference  methods 

allow  the  reconstruction  of  Inhomogeneltles  of  a  few  percent 
against  the  background  regardless  of  the  size  of  an  inhomogenelty. 
This  is  fully  supported  by  our  experimental  results  [2,24,27-293. 
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To  illustrate  the  estimation  of  the  sensitivity  of  the 
methods,  we  shall  consider  the  reconstruction  results  of  model  3, 
with  the  chain  of  irregularities  of  the  size  a  =1 00km  and 
variations  AN=0. 025*1  o’^^m"'^  (the  first  two  irregixlarities  at  the 
left)  and  AN=0.015  *10*’^"^  (the  next  three  lrr*egularities ) ,  which 
amounts  to  3-5%  of  the  maximum  Njjj=0.5*10'^^m  Due  to 
irregularities,  variations  of  the  total  pliase  are  a  few  parts  of  a 
percent  and  practically  invisible  in  fig.  19.  In  accordance  with 
the  above  estimations,  doppler  variations  are  equal  to  a  few 
percent  and  can  be  seen  easily  in  fig. 20.  The  HT  reconstruction 
results  for  model  3  are  represented  in  fig.21  (phase  RT)  end 
flg.22  (phase  -  difference  RT).  The  experiment’s  geometry  was 
taken  to  be  the  same  as  that  of  RAtE-93  (h^^l  000km,  x-j=0km, 
T2=380.4-km,  0:3=624. 91un,  T;^=809.3km).  As  an  initial  guess  we  used 
here  a  very  good  approximation  (fig. 23)  which  in  fact  coincides 
with  the  background  ionosphere.  In  spite  of  this  good  initial 
guess  the  phase  method  with  the  ART  algorithm  does  not  reveal  the 
irregularity  chain,  vdiereas  the  phase-difference  method  does 
Iridentify  the  given  structure  quite  satisfactority.  Note  that  if 
the  phase  method  is  used  with  the  MART  algorithm,  the 
iiregulai-’ities  can  be  Isolated  (fig.24).  I£4Pd’  algorithms  work 
better  within  ihe  range  of  higti  values  of  the  functions  being 
reconstructed.  However,  this  cannot  be  considered  as  an 
unconditional' indication  of  MART  algorithms  being  superior  to  ART 
ones.  If  ,\fbi* 'example,  a  similar  irregular  structure  of  even  a 
higher  intensity  ,  AN=0.04*1()*’^m"^  (fig.25)  is  located  at  a  greater 
altitude  (at '  atmt  ^6p(Dlm  not  within  the  range  of  high 

values  Of/ the  function'  being  reconstructed,  the  MRT  algorithm 
begins  distorting  strongly  the  result  (fig. 26)  in  the  attempt  to 
"attach”  such  afi  irre^ar  structure  iust  inside  the  range  of  high 


values  of  N,  in  other  words,  to  attach  the  irregularities  to  the 
maxiraum  of  the  layer.  It  goes  without  saying  that  such  a 
distortion  of  results  can  lead  to  a  wrong  interpretation  of  the 
probing  data  and  is  undoubtedly  a  serious  drawback  of  BfART 
algorithms. 

To  summarize  in  brief  the  results  of  analysing  various 
algori thins  of  solving  SLE,  it  should  be  pointed  out  that  at 
present  it  is  not  possible  to  say  which  algorithms  have 
unconditional  advantages  over  the  othes,  it  seems  to  be  impossible 
for  all  the  cases  of  ionosphere  RT.  It  is  necessary  to  establish 
the  conditions  and  areas  of  applicability  for  various  methods  and 
algori tiims  as  soon  as  possible,  therefore  now  there  is  a  vast 
field  for  examining  different  combinations  of  algorithms  as 
applied  to  problems  of  both  phase  and  pliase- difference  RT.  From 
our  point  of  view,  studies  of  rnethodn  to  solve  systems  of 
equations  [30-323  as  applied  to  the  ionospheric  RT  are  also  of 
interest.  Mansion  should  be  made  of  the  method  using 
one-dimensional  empirical  orthonormal  functions  to  reconstruct 
vertical  ionospheric  profile  [303,  the  method  based  on  the  maxiraum 
entropy  principle  [31 3 ,  a  different  orthogonal  basis  functions 
(whole  domain  functions)  reconstruction  algori tJim  [323;  several 
transfoim  techniques  [16,323.  Expansion  into  series  in  some 
continuous  basis  functions  has  both  advantages  and  disadvantages. 
In  our  opinion,  on  the  basis  of  our  several  years  of  experience  In 
experimental  RT  reconstruction,  one  is  most  unlikely  to  find  basis 
functions  to  describe  adequately  the  variety  of  ionospheric 
processes.  In  order  to  obtain  a  resolution  similar  to  that  of 
discrete  division,  the  number  of  basic  functions  imjist  be  of  the 
same  order  as  the  number  of  discrete  elements.  Besides,  the  matrix 
of  this  problem  would  be  less  rarefied  with  a  large  condition 
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number  and  the  method  Itself  is  hl^ly  sensitive  to  errors  in  the 
unknowns  coefficients  connected  with  high-order  basic  functions. 
Note  the  convenience  of  applying  the  methods  [30-32]  to  the 
phase-difference  problem  have  obvious  advantages  over  the  phase 
one.  On  the  whole  it  should  pointed  out  that  at  present  there  is 
no  one  preferable  method  to  solve  similar  tomographic  problems. 
Various  methods  may  be  better  depending  on  the  conditions  of  the 
esperlment.  It  can  be  concluded  from  this  that  it  is  necessary  to 
perform  extensive  research  to  test  vai’lous  methods  and  find  the 
optimal  methods  as  applied  to  specific  ray  RI  schemes. 

Here  we  shall  not  consider  in  detail  the  effects  of  the 
initial  quesa  on  the  reconstruction  results  and  make  some  brief 
remarlcs  only.  In  simulations  the  uniform  (in  t)  ionosphere  with 
different  levels  of  the  concentration  maxima  were  used  as  an 
initial  quess.  Changes  in  Initial  quesses  produce  small  changes 
in  the  background  level,  but  the  spatial  structure  of  the  local 
extrema  remains  the  same.  The  spatial  structure  reconstructed  by 
the  phase-difference  RT  method  can  be  said  to  weakly  depend  on  the 
initial  quess.  Generally  speaking,  the  aim  is  to  generate  an 
assembly  of  solutions  satisfying  (14)  with  a  given  accuracy 
detemined  by  experimental  errors.  Now  we  liave  developed  such 
methods  of  generatirjg  such  solution  assembly  (by  varying  the 
algorithms,  initial  quess,  etc.)  which  makes  it  possible  to  get  an 
"assembly-averaged*"  solution  and  to  estlflate  the  reconstruction 
error  distribution.  This  subject  is,  however,  too  extensive  to  be 
discussed  in  this  paper. 


3.  Program  <BOlve.for> 

This  program  solves  inverse  problem  for  Phase  RT  or 
Phase-difference  RT  by  means  of  different  algorithms:  ART,  MART, 
SIRT.  It  calculates  errors  of  reconstruction  in  metric  C  and  12  in 
dependence  of  data  errors  and  noises. 

Input  parameters  and  files: 

Parameters  <  Answer,  NP,  HR,  NRE0>  are  similar  to  saine  parameters 
of  program  <MTRIG.POR>. 

MOD.GRD  -  input  file  of  model  structure  (program  <INTEGR.POR>) 
XO.GRD  -  input  file  of  Initial  guess  (program  <INTEGR.POR> ) 

X  -  initial  guess 

Zmax  -  max  value  of  Initial  guess 

Function  <GniSS>  -  for  initial  guess 

RMAX,  RM,  Zmax,  ZSM,  B1 ,  B2  -  parameters  for  function  <GUESS> 

XIST  -  model  structure  from  input  file  <MOD.GRD> 

File  <name_F.slv>  contalnes  the  names  of  input  and  output  files: 
line  1  -  name  of  file  with  parameters  of  matric  (Fparam) 
line  2:  name1  -  name  of  input  file  with  matric 

name2  -  name  of  input  file  with  either  Phase  or  TEG  or  Doppler 
(output  files  (Pint)  from  program  <IHTEGR.FOR>) 
names  -  name  of  output  file  with  results  of  multiplications 
of  input  matric  and  reconstructed  structure 
Npi  -  array  of  contants  (2'iai)  for  each  receiver 
er_n  -  level  of  noise  in  % 

Npi,  er_n  are  Introduced  from  the  screen 

Nsolve  -  the  method  of  solution,  it  is  Introduced  from  the  screen 
REL  -  array  of  relaxation  parameter 

IMAX  -  max  number  of  nonzero  matric  elements  for  each  receiver 
AZ  -  array  of  nonzero  matric  elements  for  all  receivers 
NST  -  array  of  corresponding  column  nonzero  matric  elements  for 
all  receivers 

LST  -  array  of  number  of  all  rays  which  cross  the 

corresponding  dlscret  of  model  (SIRT  algorithm)  from 
input  file  <ST.dat>,  it  is  similar  to  LST 
in  program  <MATRIG.POR> 

Niter  -  number  of  iterations,  it  is  introduced  from  screen 
Subroutine  <DEPSyB>  -  solution  of  linear  system  equations  for  one 
ray. 

Subroutine  <ercl2>  calculates  errors  in  metric  0  and  12 


Subroutine  calculates  min  and  max  of  array 

Function  <RAN>  calculates  random  values  in  CO,  13 
Output  files: 

File  <HE0.GPD>  -  the  reconstruotion 

File  <er-3olv.dat>  -  errors  of  right  items  and  reconstruction  in 
metric  C  and  L2 

File  <Fout>  -  results  of  multiplications  of  the  matric 
and  reconstructed  structure 

Execution 
f77  solve. for 
RUN 

solve. exp 


4.  Analysis  of  the  influence  of  data  errors  and  noises  on  the 
reconstruction  results « 

Before  presenting  the  results  on  the  influence  of  errors  in 
determining  the  absolute  phase  on  the  reconstruction  results,  let 
us  estimate  the  possibilities  of  determining  the  absolute  phase. 
The  phase  method  Involves  measuring  a  linear  integral  of  the  form 
(1 )  multiplied  by  constant  of  the  order  of  unity  E23,  which  is 
insignificant  here.  The  basis  for  difficulties  appearing  In  the 
determination  of  the  linear  integral  (6)  is  that  the  phase  value 
is  very  high.  For  typical  maximum  values  N  "  10^^  m”^,  X=2m  and  a 
ray  length  in  the  ionosphere  of  the  order  of  a  thousand 
kilometers,  (p  Is  thousands  of  radian.  Thus,  the  problem  arises  of 
isolating  the  "Initial  phase"  ({>q  =2^1,  which  must  be  added  to  the 
measured  (within  Z%)  A^)  to  obtain  the  complete  phase  ({>  =  c{)q  +  Acj) 
or  the  linear  integral  (6).  To  explain  the  difficulties  arising, 
we  will  examine  the  possibility  for  isolating  the  initial  phase  in 
the  presence  of  minor  horizontal  gradients. 

We  represent  the  concentration  in  the  form  of  an  expansion, 
where  the  regular  spherlcally-syrametric  background  N^Cr)  is 
Isolated  N(r,a)  =  NQ(r)  +  N’ (a  -  cx^),  N’ (r)  s  ;  a^((j))  ia 


3^ 


the  angle  of  Intersection  of  the  ionospheric  masiimm  by  the  ray; 
its  Ticinity  malies  the  primary  contribution  to  integral  (6).  In 
this  case,  retaining  the  first  teims  of  the  expansion  in  terms  of 
powers  h/R  in  (6),  we  obtain 


XT, 


COS(|) 


tg^-(|) 


h 


J  M  (h)(lh  -  Xr„  - J  -  No(h)dh  + 

®  COSCj)  ^ 


cosv!) 


J  N'(h)(a(h)-0(^)dh 


(22) 


If  there  were  no  horizontal  gradient,  it  might  be  possible  to 
measure  A(})(tj))  for  various  angles  <]>  in  the  range  h<p  and  obtain  a 
system  of  linear  equations  from  which  it  would  be  possible  to 


determine  j  N^dh,  J  dh  h  N^/H 
functional  dependence  on  (|)  would 
TEC  and  other  momenta  of  the  funo 


etc.  In  other  words,  the  known 
malrce  it  possible  to  isolate  the 
tlon  N^Ch).  nowevef,  the  presence 


of  the  teiTii  with  N’  greatly  complicates  the  situation  when  its 


value  becomes  comparable  to  2%.  We  will  estimate  this  term  in  the 
case  of  nearly  vertical  sounding  |al,  ja^t  1;  for  this,  in  the 
integrand  we  replace  a(h)-  by  4)(h-hi^^)/(R+h)  4)(h-hjy^)/R  ,  h^ 
is  the  height  of  the  maximum;  this  asymptotic  equation  follows  (at 
low  angles)  from  (4):  (a--a^)/(1~R/r)  =  (a-a^^)  (R+h)/h. 
Therefore,  such  methods  for  isolating  the  constant  component  will 
be  suitable  under  the  condition 


N’(h)(h-li)  , 

A(|)  r  - Tn_  dh  «  29C. 

®  R 


(23) 


The  typical  values  of  dlVda  in  the  presence  of  a  "trough”  in 
the  ionosphere  dWda  10*^^  radian"**,  then  (23)  is  only 
satified  for  Atj)  «  10“^.  However,  at  angles  of  fractions  of  a 
degree  it  is  practically  impossible  to  determine  the  functional 
dependence  on  (j)  in  the  presence  of  noises.  Inequality  (23)  applies 


to  the  case  of  deteimination  of  the  TEC  along  the  Yertical.  The 
llmltatloPx  on  the  horizontal  gradient  hecomes  even  more  strict  In 
oblique  soimding.  A  detailed  analysis  of  a  number  of  traditional 
techniques  for  determination  the  constant  phase  was  proYlded  In 
[333.  These  techniques  Include  combined  tecimiques  Doppler  and 
Faraday  measurements  simultaneously.  The  results  of  numerical 
modeling  [333  showed  tlmt  in  the  presence  of  a  trough  and 
characteristic  gradients  dWda  10'‘^  m“^  radian"'’  the  error  In 
determination  of  the  constants  Yarles  within  100-1000%!  This 
agrees  with  estimate  (23)  and  Indicates  that  practical 
deteiminatlon  of  the  linear  integral  of  the  electron  concentration 


is  unrealistic  In  the  pr-asenGS 


of  characteristic  horizontal 


gradients  iii  the  ionoaphere,  while  it  is  specifically  this  case 


which  is  of  interest  for  RT. 


There  is  one  more  method  for  determining  the  TEG  in  the 


presence  of  horizontal  gradients  [34-363  which  is  baaed  on 
recording  of  satellite  signals  at  a  pair  of  separated  recelYlng 
stations.  With  data  on  the  phases  on  a  fixed  base,  it  is  possible 
to  compare  the  pair  of  linear  equations  with  the  nearly  identical 
last  terms  of  (22),  l.e.,  the  rays  traveling  to  different 

receivers  should  intersect  the  ionospheric  maxiinum  at  one  point. 
L13cewlse,  a  pair  of  equations  for  another  moment  In  time  leads  to 
a  system  for  the  initial  phase,  in  which  the  Influence  of  the 
gradient  term  will  be  reduced.  Otherwise,  condition  (23)  is 
replaced  by  a  less  strict  condition.  However,  it  is  impossible  to 
completely  eliminate  the  Influence  of  horizontal  gradient  and 
subsequent  derivatives  on  the  result  of  determination  of  the 
initial  phase  by  a  such  a  method.  The  methods  used  [28-303  have  an 
error  by  constant  of  about  10%  or  at  least  a  few  percents.  Of 
coiiTse,  it  is  always  possible  to  propose  a  variant  of  the 


recording  method  in  which  it  is  possible  to  determine  the  initial 
phase;  but  to  do  so,  one  must  perform  a  multifrequency  reception 

using  an  array  of  seTeral  reoeiTers, 

Thus,  by  virtue  of  the  natiu^e  of  the  phase  measurements,  it 
is  inadvisable  to  reduce  the  problem  of  ionospheric  RT  to  a 
problem  in  linear  integrals.  Determination  of  initial  phase  by  the 
simplest  recording  systems  leads  to  major  errors,  ?Mle  the  use  of 
complex  fflultiposltion  and  nroltifrequency  systems  is  not  Justified 
here,  since  a  different  solution  of  the  problem  is  possible  by  the 


phase-difference  or  Doppler  measurements  without  determination  of 
the  initial  phase. 

Here  are  the  results  of  numerical  simulation  of  the 
reconstruction  of  the  Ionosphere  section  for  tj/pical  errors  in 
finding  the  absolute  phase  (TEC),  A  simple  model  of  the  ionosphere 
with  a  trough  and  a  positive  inTiomogeneity  at  the  left  edge  of  the 
trough  (model  1 )  is  used  for  nrmerlcal  simulation.  It  was  assumed 


that  satellite  radio  probing  (hQ=1000  Ion)  was  performed  at  the 
frequencies  400  and  150  MHs  (A.=£ffi)  snd  the  receivers  are  located 
at  the  points  with  the  coordinates  a\j=0  lurij  t:2==4T5  km,  i^=14351an. 
A  homogeneous  ionosphere  having  no  trough  was  used  as  a  initial 


guess  Pig.27.  Fig.  28  and  Fig.  29  show  the  results  of  the 
reconstruction  with  the  help  of  the  phase  RT  with  -3%  and  -5% 


accuracy  in  deteming  the  TEG  (which  corresponds  the  error  In  the 
value  of  the  absolute  of  6%  and  lOrc),  the  signs  of  the  errors 
don’t  coincide  for  different  receivers.  These  figures  illustrate 


an  extremely  poor  quality  of  reconstruction  using  the  phase  RT 
method  with  typical  errors  in  determining  the  TEC:  even  the 
principal  features  of  this  simple  model  structure  are  not 
recovered  and  at  the  same  time  some  heav^*'  artefacts  are  present. 
Fig.30  shows  the  dependences  of  the  reconstruction  errors  6^ 


and  6^  on  the  error  In  determining  the  absolute  phase  (-21™), 
on  the  number  m.  It  can  be  seen  that  even  errors  of  a  few  units  2% 
lead  to  a  poor  quality  of  reconstruction.  It  should  be  noted  that 
if  the  SIRT  algorithm  is  used,  as  it  was  done  by  other  authors  in 
[12,17,19],  rather  than  the  ART  algorithm,  the  reconstruction 
results  similar  to  those  in  Figs. 28-29  will  at  first  sight  be 
better.  They  haye  a  more  regular,  less  "chaotic”  character.  The 
averaging  SIRT  algorithms  result  in  the  reconstruction  being 
weaker  dependent  on  measiirement  errors,  i.e.'  the  SIRT  algorithms 
are  less  sensitive  to  small  variations  of  data  and  therefore 
poorly  reflect  the  fine  structure  of  the  reconstructed 
cross-sections.  The  SIRT  algorithms  can  reconstruct  the  general 
background  quite  satisfactorily  but  often  they  fall  to  reveal  even 
the  existing  trough  if  it  is  20-40%  smaller  than  the  maximum.  If 
the  signs  of  the  errors  in  determining  the  absolute  phase  are  the 
same  for  all  the  receivers,  the  reconstruction  quality  is  somewhat 
better.  Nevertheless  at  the  level  of  errors  greater  then  10%,  the 
reconstruction  quality  is  still  poor.  In  Fig.31  we  represent  the 
result  of  reconstruction  of  the  model  1  using  the  phase- 
difference  RT  method  (the  matrix  was  built  by  piecewise-planar 
approximation).  One  cmi  easily  see  that  the  principal  features  of 
the  ionosphere  section  are  reconstructed  quite  well.  Numerical 
simulation  of  the  reconstruction  of  various  ionospheric  structures 
carried  out  by  us  proved  a  noticeable  advantage  of  the 
phase-difference  RT  method  over  the  phase  RT  with  typical  errors 
in  determining  the  absolute  phase. 

Note  another  significant  limitation  of  RT  methods  connected 
with  deviation  of  receivers  from  the  satellite  flight  plane. 
Similar  deviation  can  also  lead  to  sigiificant  errors  in 
determining  the  absolute  phase.  Let  iis  Introduce  the  distance  p  in 


the  direction  perpendicular  to  this  plane,  and  the  deviation  angle 
tg7=p/L  where  L  is  the  (tilted)  distance  to  the  satellite.  Then  it 
can  easily  be  shovm  that  the  difference  between  the  pliase  detected 
by  the  receiver  in  the  flight  plane  and  that  detected  by  receiver 
off  the  plane  by  p  will  be: 

A0  «  ?irgN^(1/cos7  -  1 )+  ;\,rgSln7  <aN/ap>  L^/2  (24) 

L 

2  r  aN 

where  <aN/ap>  =  — ^  -  Idl  -■  the  ray  "averaged" 

^  i  p=o 

transverse  gradient  of  the  electron  density,  N^.  is  the  TEG  along 
the  ray.  The  correction  associated  with  changing  tilted  distance 
(the  first  suinmand  in  (24))  cen  easily  be  taken  into  accomt. 
However,  strong  transverse  gradients  can  noticeably  change  the 
phase,  and  for  the  phase  method  the  following  condition  must  be 
satisfied: 

K  r^  sln7  <“|p->  «  2'/t  (25) 

with  the  typjlcal  gradient  <aN/ap>~  10^  and  distances  L"'10m^ 
this  fact  results  In  a  strict  3  Imitation  of  p«10km.  However 
this  limitation  was  not  taken  into  account  in  the  experiments 
[17,19,203.  The  Bheme  of  experiment  Cl 7 3  seems  to  be  strange  and 
surprising  because  the  Kallntngrad-RlgB-Lenlngrad  track  makes  an 
angle  of  the  order  of  45®  with  the  satellite  trajectory 
projection.  The  angle  between  Kinma-Oulu  direction  [193  and  the 
satellite  trajectory/  projection  Is  also  greater  then  40®.  Note 
that  the  deviation  of  the  receiver  from  the  satellite  plane  only 
due  to  the  Earth's  rotation  during  the  recording  time  of  the  order 
of  10  minutes  can  reach  a  hundred  kilometers  in  the  middle 
latitudes,  which  is  quite  significant  for  phase  RT.  For  the 


phase-difference  RT  it  is  only  relative  smallness  of  the  change  A(}) 
in  respect  to  the  basic  phase  that  is  required,  which  leads  to  the 
inequality  <5K^/0p>p  «  tliat  is  practically  always  fulfilled. 

Consider  the  Influence  of  the  noise-type  errors  on  the 
reconstraction  results.  The  physical  origin  of  such  noises  may  be 
connected  both  with  instrumental  errors  and  external  noises. 
Methods  of  phase-difference  ET  provide  quite  satisfactory  results 
even  with  considerable  errors  in  experimental  data.  Pig.32  shows 
the  reconstruction  by  the  method  of  phase-difference  RT,  but  in 
this  case  the  data  with  the  noise  level  being  10%  of  maximum 
amplitude  of  Doppler  data  were  used.  One  can  see  that  random 
measurements  errors  slightly  affect  the  reconstruction  results, 
almost  there  is  no  defference  between  Fig. 31  and  Fig. 32.  This 
conclusion  is  confirmed  by  fig. 33  showing  the  dependence  of  the 
reconstruction  error  in  the  1^  metrics  on  the  noise  level.  Even 
the  50%  noise  level  changes  the  reconstruction  error  only 
slightly,  which  can  be  explained  by  mutual  compensation  and 
effective  "averaging"  of  noises  in  the  process  of  tomographic 
reconstruction.  In  the  phsae-difference  RT  experiments 
£24,27-291  the  noise  level  does  not  exceed  a  few  percent.  The 
results  of  the  simulation  also  show  that  the  influence  of  noises 
on  the  phase  RT  method  is  rather  weak,  because  within  the  frame 
work  of  this  model  it  Is  possible  to  consider  weak  noises  within 
the  27c-lnntemal  only,  which  amounts  to  10~®  in  typical 
ionospheric  conditions.  Thus,  for  the  phase  RT  method,  errors  in 
determining  the  "initial  phase"  are  of  paramount  Importance  and  it 
is  difficult  to  avoid  such  errors  due  to  numerous  factors,  such  as 
horisontal  gradients,  deviations  of  receivers  from  the  plane  of 
recording,  changes  in  the  position  of  receivers  owing  the  Earth’s 
rotation,  etc. 
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Table  1 

Direct  problem  errors 


Model  5  Model  2  Model  7 


Operator  type 

Sco 

«2 

82 

A 

0.0344 

0.0087 

0.1241 

0.0803 

0.3160 

0.2105- 

B 

0.01 4T 

0.0065 

0.0933 

0,0540 

0.0500 

0.0323 

0 

0.0101 

0.0071 

0.0898 

0.0534 

0.0468 

0.0303 

D 

0.0005 

0.0002 

0.0449 

0.0244 

Table  2 

Inverse  problem  errors 


Model 

5 

Model  2 

Model  7 

OpQ|»atpj 

?  trpe  5^ 

02 

5 

CO 

: 

A 

0.26 

0,21 

0.328 

0.225 

Q.407 

0.327 

B 

D.1T 

0,15 

0.317 

4.215  ^ 

0,180 

0.127, 

0 

0.16 

QrU 

0.316 

0.211 

0.118 

h 

0.16 

0.16 

0,080 

0.073 

Table  3 

Inverse  problem  errors 


Operator  type 

Model 

6 

^2 

Model  2 
^0  '  ^2 

Model  7 

^00  ^2 

A 

0.25  ? 

0.22 

0.421 

0406 

0,660 

0.362 

B 

0.16  ? 

0,14 

0.339 

0t#6 

p.jto 

0.089 

C 

0.15  7 

0.14 

0.339 

04^4 

lifti 

0.075 

D 

0.16  7 

5.13 

0.338 

0.224 

Orili 

0.069 

